Одна зависимая (объясняемая) переменная: y
Несколько регрессоров (предикторов, объясняющих переменных): x, z…
По каждой переменной n наблюдений: \(y_1, y_2...y_n\)
Модель — это формуля для объясняемой переменной.
Возьмём например данные по машинам 1920-х годов. Тут видимая линейная взаимосвязь.
Модель для этих данных может иметь вид \(y_i = \beta_1 + \beta_2x_i + \varepsilon_i\)
Что здесь есть что? У нас есть наблюдаемые переменные,
y — это длина тормозного пути
x — это скорость, с которой ехала машина.
Есть неизвестные коэффициенты $ _1, _2$
случайная составляющая, ошибка
То есть, \(\beta_2\) показывает — насколько увеличивается тормозной путь, если машина разгонится на один лишний километр в час. И есть некая случайная составляющая \(ε_i\), это может быть все что угодно:
то есть это та часть, по которой у нас нет возможности предсказать, но, тем не менее, вот эта случайная ошибка она входит в y.
Придумать адекватную модель
Получить оценки неизвестных параметров \(\widehat\beta_1, \widehat\beta_2\)
Спрогнозировать, заменив неизвестные параметры на оценки \(\widehat y_i = \widehat\beta_1 + \widehat\beta_2 x_i\)
Как найти \(\beta_1, \beta_2\)? Собственно методом наименьших квадратов.
Если мы придумали какие-то оценки, \(\beta_1, \beta_2\) то соответственно возникает такое понятие как ошибка прогноза
\[\widehat \varepsilon_i = y_i - \widehat y_i\]
Есть суммарная ошибка, чтобы суммарная ошибка не занулялась одна в плюс, другая в минус, не компенсировали друг друга, мы возведем в квадрат. И посчитаем сумму квадратов ошибок прогноза, то есть сумма \(\widehat \varepsilon_i ^2\)
\[Q(\widehat\beta_1, \widehat\beta_2) = \sum_{i=1}^{n} \widehat\varepsilon^{2}_{i} = \sum_{i=1}^{n} (y_i - \widehat y_i)^2 \]
Суть метода МНК
Возьмите в качестве оценок такие коэффициенты \(\widehat\beta_1, \widehat\beta_2\) при которых сумма квадратов ошибок прогноза будет минимальна
Всё сводится к тому, что в модели \[M_1: y_i = \beta + \varepsilon_i\]
\[\widehat\beta = \bar y\]
Случай для одного регрессора
Тут несколько сложнее, подробности на скриншоте
Случай для двух регрессоров
Что мы получили:
Случай для двух регрессоров
Ещё раз пробежимся по терминам
Случай для двух регрессоров продолжение
Изобразим на графике, всё то, о чем мы только что проговорили.
Графический вывод
В частности \(\widehat\beta_1\), т.е. точка на оси ординат от пересечении с прямой регрессии, называют пересечением или Intercept-ом
Метод наименьших квадратов подбирает прямую так, чтобы суммарные расстояния от точек до прямой были минимальными.
Графический вывод продолжение
Случай множественных регрессоров принципиально не отличается от двух регрессоров. Поэтому рассмотрим на примере трёх предикторов.
Множественные регрессоры
Суммы квадратов
Что есть что:
RSS (Residual Sum of Squares) — этот показатель меряет, насколько велики эпсилон с крышкой, насколько они далеко лежат от нуля.
TSS (Total Sum of Squares) — этот показатель меряет, насколько каждый из \(y_i\) не похож на \(y\) среднее. Если \(y_i\) далеко лежит от \(y\) среднее, соответственно, это слагаемое в сумме квадратов будет большим. И вся сумма квадратов будет большой.
ESS (Explained Sum of Squares) — она показывает, насколько прогнозное значение \(y_i\) с крышкой далеко легло от \(y\) среднего.
Язык эконометрики во многом — это язык линейной алгебры, нужно его знать.
Обозначения
Маленькой буквой \(y\) мы будем обозначать вектор, то есть столбик из всех игреков, записанных друг под другом: у_1, у_2 и так далее, у_n.
Ну, соответственно, \(х\) маленькое — это все иксы: х_1, х_2 и так далее, х_n, записанные друг под другом.
Аналогично \(\widehat\beta\).
И еще введем обозначение: единичку со стрелочкой. Это будет вектор-столбец, то есть столбик из единичек в количестве n штук, потому что у нас в модели будет n наблюдений.
Тогда для нашей модели
\[\widehat y = \widehat\beta_1 \cdot \overrightarrow 1 + \widehat\beta_2 \cdot x + \widehat\beta_3 \cdot z\]
Таблица регрессоров
Длина вектора
Где у нас бывают длины векторов:
Пример длины вектора
Есть одно замечательное применение. Если скалярное произведение векторов равно нулю, значит они перпендикулярны.
Скалярное произведение векторов
x <- c(1,2, -3)
y <- c(-6,0,-2)
sum(x * y)
## [1] 0
n-мерное простарнство
Есть вектор y, есть вектор из одних 1, я продолжаю этот вектор до прямой и оказывается, что
Соответственно, мы получили попутно еще один факт: если любой вектор проецировать на вектор из 1, то получится вектор средних значений
Напомню, что мы вывели условие первого порядка
Геометрическая интерпретация условий первого порядка
Это означает, что вектора перпендикулярны.
Облачко это все те вектора, которые можно получить, складывая с некоторыми весами вектор x, вектор z и вектор из единичек. Это гиперплоскость.
Мы получаем следующую интрпретацию метода наименьших квадратов:
Первый геометрический факт Прогнозы, которые мы получаем по методу наименьших квадратов — это проекция исходного вектора зависимых переменных y на множество векторов, получаемых с помощью сложения с разными весами вектора из единичек, вектора x и вектора z.
Второй геометрический факт Если я спроецировать вектор y на прямую, порождаемую вектором из единичек, и спроецировать y с крышечкой на ту же самую прямую, то эти проекции попадут в одну и ту же точку.
\(TSS = ESS + RSS\)
\(\frac{ESS}{TSS}= \frac{BC^2}{AB^2}=(\frac{BC}{AB})^2= (cos \varphi)^2\)
Геометрчиеские факты
Если в регрессию включён свободный член \(\beta_1\) то действуют следующие правила:
Правила при включении свободного члена в регрессию
Эти правила позволяют придумать простой показатель качества работы модели — коэффициент детерминации.
Чем прогнозы точнее похожи на настоящие значения, тем меньше будут ошибки прогнозов и тем меньше будет сумма квадратов ошибок прогнозов RSS. Соответственно \(\frac{ESS}{TSS} = R^2\) будет примерно равен 1 если RSS будет у ноля. Соответственно мы получим коэффициент детерминации, который всегда лежит от 0 до 1.
Другими словами — коэффициент детерминации это доля объяснённого разброса в общем разбросе.
Правила при включении свободного члена в регрессию
С одной стороны, коэффициент детерминации — это доля объясненной дисперсии игрек, доля объясненного разброса.
С другой стороны, коэффициент детерминации — это выборочная корреляция между прогнозами и настоящим игрек, взятое в квадрат.
Чем коэффициент детерминации выше, тем больше предсказание похожи на реальные значения. Чем коэффициент детерминации выше, тем выше доля объясненной дисперсии.
Посмотрим как зависит фертильность от других
t <- swiss
# нарисовать много диаграм рассеивания
# install.packages("GGally") # матрица диаграмм рассеяния
library("GGally")
## Loading required package: ggplot2
str(t)
## 'data.frame': 47 obs. of 6 variables:
## $ Fertility : num 80.2 83.1 92.5 85.8 76.9 76.1 83.8 92.4 82.4 82.9 ...
## $ Agriculture : num 17 45.1 39.7 36.5 43.5 35.3 70.2 67.8 53.3 45.2 ...
## $ Examination : int 15 6 5 12 17 9 16 14 12 16 ...
## $ Education : int 12 9 5 7 15 7 7 8 7 13 ...
## $ Catholic : num 9.96 84.84 93.4 33.77 5.16 ...
## $ Infant.Mortality: num 22.2 22.2 20.2 20.3 20.6 26.6 23.6 24.9 21 24.4 ...
ggpairs(t)
Уже тут можно посмотреть много всяких корреляций.
Перейдём к оценке
model2 <- lm(data = t, Fertility ~ Agriculture + Education + Catholic)
# КОэффициенты
coef(model2)
## (Intercept) Agriculture Education Catholic
## 86.2250198 -0.2030377 -1.0721468 0.1452013
# Спрогнозированные значения
fitted(model2)
## Courtelary Delemont Franches-Mnt Moutier Neuveville
## 71.35382 79.73758 86.36550 76.21257 62.05992
## Porrentruy Broye Glane Gruyere Sarine
## 84.70365 77.94869 77.98965 82.07990 76.37831
## Veveyse Aigle Aubonne Avenches Cossonay
## 81.01451 62.00804 65.34456 61.67811 67.20324
## Echallens Grandson Lausanne La Vallee Lavaux
## 72.85406 71.22373 54.02437 62.00809 62.16632
## Morges Moudon Nyone Orbe Oron
## 64.12130 72.47751 65.22299 69.41765 71.04507
## Payerne Paysd'enhaut Rolle Vevey Yverdon
## 66.61076 70.48740 64.27982 63.09324 68.48321
## Conthey Entremont Herens Martigwy Monthey
## 81.11782 77.02791 80.38838 78.28372 84.09311
## St Maurice Sierre Sion Boudry La Chauxdfnd
## 75.54878 80.27332 73.53528 66.37864 74.87034
## Le Locle Neuchatel Val de Ruz ValdeTravers V. De Geneve
## 70.52554 50.79966 71.80743 76.17918 35.30542
## Rive Droite Rive Gauche
## 52.99371 57.97821
# Остатки
residuals(model2)
## Courtelary Delemont Franches-Mnt Moutier Neuveville
## 8.8461774 3.3624191 6.1345049 9.5874340 14.8400828
## Porrentruy Broye Glane Gruyere Sarine
## -8.6036475 5.8513084 14.4103471 0.3201012 6.5216935
## Veveyse Aigle Aubonne Avenches Cossonay
## 6.0854871 2.0919628 1.5554442 7.2218873 -5.5032423
## Echallens Grandson Lausanne La Vallee Lavaux
## -4.5540632 0.4762715 1.6756342 -7.7080933 2.9336804
## Morges Moudon Nyone Orbe Oron
## 1.3786987 -7.4775133 -8.6229883 -12.0176461 1.4549265
## Payerne Paysd'enhaut Rolle Vevey Yverdon
## 7.5892409 1.5125978 -3.7798150 -4.7932370 -3.0832083
## Conthey Entremont Herens Martigwy Monthey
## -5.6178154 -7.7279097 -3.0883806 -7.7837172 -4.6931098
## St Maurice Sierre Sion Boudry La Chauxdfnd
## -10.5487835 11.9266828 5.7647206 4.0213575 -9.1703410
## Le Locle Neuchatel Val de Ruz ValdeTravers V. De Geneve
## 2.1744592 13.6003353 5.7925740 -8.5791790 -0.3054172
## Rive Droite Rive Gauche
## -8.2937095 -15.1782122
# ПОказатель RSS
deviance(model2)
## [1] 2567.884
# Показатель R квадрат
report <- summary(model2)
report$r.squared
## [1] 0.6422541
# Это тоже самое что
cor(t$Fertility, fitted(model2))^2
## [1] 0.6422541
Вторая неделя — стат.свойства оценок коэффициентов
Формулируем стандартные предпосылки — они ведут к свойствам — свойства позволят строить доверительные интервалы.
Для формулировки предпосылок — сформулируем поняти условное мат.ожидание.
Условия мат. ожидания
Если формально: это некая случайная величина s с тильдой, которая является, с одной стороны, функцией от r, и, с другой стороны, эта самая s с тильдой очень похожа на s, а именно: математической ожидание от s с тильдой равно математическому ожиданию от s, и ковариация между s и любой функцией от r равна ковариации между s с тильдой и той же самой функцией от r.
То есть получается, что s с тильдой и s очень похожи, s с тильдой и s невозможно отличить, если смотреть только на математическое ожидание или на ковариацию с r, с r-квадрат, с любой функцией от r
Условное мат. ожидание — это Верно! \(E(y_i | x_i)\) это наилучший прогноз \(y_i\) формулируемый с помощью \(x_i\).
На практике
Условия мат. ожидания
Здесь показано как рассчитывать
Условия мат. ожидания
Если величины непрерывны и есть совместная функция плотности
Условия мат. ожидания
Условия мат. ожидания
Найдите \(E(x_i^2+2x_i | x_i)\) это \(x^{2}_i+2x_i\). Потому что если мы знаем \(x_i\), то и всё выражение слева мы легко спрогнозируем.
Условия мат. ожидания
Свойства:
Условия мат. ожидания
Пример
Упростите \(Var(x^{2}_i+2x_i+\varepsilon_i|x_i)\) это равно \(Var(\epsilon_i | x_i)\) Потому что Вспомним свойство \(Var(s+h(r) | r)=Var(s|r)\) В роли s выступает \(\varepsilon_i\). А если интуитивно: хотим спрогнозировать выражение слева, зная \(x_i\) часть \(x_i^2+2x_i\) нам доподлинно известна и на точность прогнозирования (а именно её меряет условная дисперсия) не влияет.
-
Запомнить
Интерпретируем дисперсию как квадрат длины случайной величины,
Интерпретируем корреляцию между случайными величинами как косинус угла между ними.
С помощью условного математического ожидания мы сформулируем стандартные предпосылки на случайную составляющую \(\varepsilon\), а миеннто три предпосылки
Математическое ожидание от каждой случайной составляющей при известных иксах, то есть при известном каждом регресссоре для каждого наблюдения, я пишу коротко матрицу X. Условное математическое ожидание от \(\varepsilon_i\) при всех известных регрессорах равна 0
Условное математическое ожидание от \(\varepsilon^2\) при всех условных регрессорах равна \(\sigma^2\). Или можно точно так же сказать, что условная дисперсия \(\varepsilon_i\) при известной матрице X равна \(\sigma^2\).
Kовариация между \(\varepsilon_i\) и \(\varepsilon__j\) при фиксированной матрице X равна нулю.
Предпосылки
Предпосылки про ковариацию и дисперсию можно коротко записать с помощью такого понятия как ковариационная матрица
Предпосылки
Когда мы говорим «ковариационная матрица некоего случайного вектора», мы имеем в виду здоровую табличку чисел. Первое число в первой строке — это дисперсия \(\varepsilon_1\) . Второе число в первой строке, то есть первая строчка, второй столбец, (1,2) координаты, — это ковариация \(\varepsilon_1\) и \(\varepsilon_2\).
Соответственно, в ковариационной матрице, скажем, в третьей строчке, втором столбце находится ковариация \(\varepsilon_3\) и \(\varepsilon_2\), это третья строчка, второй столбец. С
оответственно, в этой матрице находятся и все дисперсии каждого \(\varepsilon\) , и все попарные ковариации: \(\varepsilon_1\) с \(\varepsilon_3\) , \(\varepsilon_2\) с \(\varepsilon_7\) … все-все-все ковариации и дисперсии находятся в одной матрице, в одной табличке чисел.
Соответственно, первые наши две предпосылки на ε можно как сформулировать? Что ковариации равны нулю, а на диагонали дисперсии равны \(\sigma^2\). То есть матрица принимает такой простой вид.
Предпосылки
Соответственно, в нашем случае наши предпосылки можно записать как ковариационная матрица вектора \(\varepsilon\) при фиксированных регрессорах X равна \(\sigma^2\) умножить на эту самую единичную матрицу, которая обозначается буковкой I, сокращение от английского «identity».
У нас есть предпосылки, что дисперсия \(\varepsilon\) при фиксированных регрессорах X равна \(\sigma^2\) умножить на единичную матрицу, что на самом деле просто означает, что дисперсия \(\varepsilon_i\) при фиксированном X равна \(\sigma^2\), а ковариация разных \(\varepsilon\) при фиксированном X равна нулю.
И у нас есть линейная модель,
Итоги
Это для примера двух объясняющих переменных. И, соответственно, эти предпосылки позволяют посчитать дисперсию любой оценки МНК \(\beta_j\) с крышкой и любую ковариацию \(\beta_j\) с крышкой и \(\beta_l\) с крышкой.
Мы посчитаем для начала дисперсию, условную, оценки метода наименьших квадратов и ковариацию оценки метода наименьших квадратов для случая парной регрессии.
Итоги
Условия нашей задачи
Что мы получаем. Запишем примечания
Примечания
Посчитаем в рамках предположений ковариацию \(\widehat\beta_2\) и среднее
ПОлучаем не что иное как
Примечания
Итого в парной регрессии мы имеем
Итого в парной регрессии
Если регрессор \(z\) сильно коррелирован с другими регрессорами, то Величина \(RSS_z\) будет примерно равна 0 и поэтому дисперсия \(Var(\widehat\beta_z|X)\) будет большой
Теорема
Штрих — это транспонированная матрица.
Свойства
Этот фрагмент особо полезен тем, кто знает линейную алгебру. Оказывается, средствами линейной алгебры легко не только сразу посчитать дисперсию \(\widehat\beta_1\) с крышкой или \(\widehat\beta_2\) с крышкой, одним махом можно найти все эти дисперсии и все ковариации.
Свойства ковариационных матриц
Свойства
В этой матрице находятся и дисперсии условные каждого \(\widehat\beta\), и ковариации каждого \(\widehat\beta\) с каждым \(\widehat\beta\). И в скалярном виде мы их выводили долго по отдельности и только для случая парной регрессии, а в матричном виде они выводятся легко, все и сразу, и этой формулой можно будет пользоваться.
Свойства
Константа \(\sigma^2\) неизвестна, и именно эта константа входит в формулу для условной дисперсии любого оценённого коэффициента. А нам хочется строить доверительные интервалы для коэффициентов, проверять гипотезы, поэтому нам нужен какой-то способ оценить неизвестную константу.
Свойства
Оценки коэффициентов, которые мы получаем методом наименьших квадратов, обладают рядом замечательных статистических свойств. Эти свойства, их очень много, и поэтому, чтобы как-то их было проще все осознать, мы их поделим на три части
Свойства
Предпосылки для применения метода наименьших квадратов и для получения хороших свойств коэффициентов.
Предпосылки
Свойства
Предположения на \(\varepsilon_i\)
Свойства
Предпоссылки на регрессоры
Свойства
Когда все эти предпосылки, которые мы сформулировали:
когда все они выполнены, мы получаем, что верны три группы свойств.
Свойства
Это очень хорошее свойство. Это говорит о том, что наш метод, конечно, может ошибаться, и оценка \(\widehat\beta\), которую мы получаем, может не совпадать с настоящим \(\beta\), но \(\widehat\beta\) иногда будет больше настоящего \(\beta\), иногда будет меньше настоящего \(\beta\), но в среднем,\(\widehat\beta\) попадает то влево, то вправо, но в среднем попадает в неизвестный коэффициент \(\beta\)
Свойства
Это очень хорошее свойство. Оно говорит о том, что если вы хотите простую оценку, то есть линейную, хотите оценку несмещенную, которая в среднем бы попадала в неизвестную истину, то ничего лучше оценок метода наименьших квадратов у нас не получится. То есть математически это означает, что условная дисперсия\(\widehat\beta\) при фиксированных регрессорах альтернативного больше либо равна, чем условная дисперсия \(\widehat\beta\), полученная по методу наименьших квадратов опять же при фиксированных регрессорах.
Свойства
То есть это свойство говорит, что оценка \(\widehat\sigma^2\), предложенная нами, — RSS, делённая на (n – k), — она тоже несмещённая, и она несмещённым образом оценивает \(\sigma^2\).
Это те, что предполагают что число экспериментов велико
Свойства
Свойства
Проверять гипотезы можно в двух случаях
Число наблюдений велико
Случайные ошибки нормальны
Свойства
С ростом стандартной ошибки ширина доверительного интервала для коэффициента увеличивается
Предпосылки теста
\(H_0\) — проверяемая гипотеза
\(H_1\) или \(H_a\) — альтернативная или другими словами конкурирующая гипотеза
Формула для вычисления статистики
Закон распределения статистики при верной \(H_0\)
Формулируем гипотезу \(H_0\) и выбираем уровень значимости. Это вероятность ошибки первого рода или вероятность отвергнуть \(H_0\) при условии, что она верна. \(\alpha = P (отвергнуть H_0| принять H_0)\). Можно формулировать гипотезу и под неё собирать данные
Рассчитываем наблюдаемое значение тестовой статистики \(S_obs\)
Находим критическое значение тестовой статистики \(S_cr\)
считается это всё вот так
Пример
qt() — это функция в R
Мы опираемся на теорему RSS делённое на сигма квадрат
qchisq() — это функция в R
Пример
Предпологаем что коэффициент при доле сельскохозяйственного мужского населения равен нулю — другими словами фертильность не зависит от показателя того, насколько этот регион является сельскохозяйственным. Есть три способа проверить эту гипотезу
Пример
В литературе, скажем в научных статьях, очень часто стандартные ошибки выписывают под коэффициентами.
При оценке линейной модели регрессии, любой статистический пакет выдает более-менее стандартную табличку. Вот такую табличку выдает R.
Пример
Первый столбик в ней это, собственно, оценки коэффициентов. Т.е. смысл этого столбика, что мы можем записать уравнение линейной регрессии по используя эти коэффициенты.
Второй столбик — это стандартные ошибки. Это корни из диагональных элементов ковариационной матрицы
Третий столбик — это компьютер автоматом проверяет гипотезу о том, что на самом деле, зависимости от данной переменной нет. Он делает это при помощи Т-статистики. Т.е. это первый столбец делить на второй. T-статистика, проверяющая гипотезу о незначимости коэффициента может принимать любое значение. Знак T-статистики определяется знаком оценки коэффициента, и по модулю она может быть произвольной
P-value
Если \(H_0\) не отвергается, это говорит о том, что зависимости нет.
Надо говорить аккуратно: «\(H_0\) не отвергается» — это означает, что данные не противоречат гипотезе \(H_0\)
Значимость и существенность. Значимость это статистический факт — равен коэффициент нулю или не равен. Существенность — это на сколько он не равен нулю.
Стандартизировать коэффициенты перед анализом. Чтобы получить одинаково интерпритируемые единицы измерения
К сожалению, очень часто распространена такая порочная практика, что исследователь берет, включает кучу, не задумываясь о теоретической модели, включает кучу объясняющих переменных в свою модель и выбирает те из них, которые по t-статистикам оказались значимы. Это подход неправильный, поскольку как только мы согласились на некую вероятность ошибки первого рода, например, мы выбрали вероятность ошибки первого рода типичную в экономических приложениях 5%. «Я запустил регрессию на кучу переменных и отобрал те, которые значимы» — это неправильный подход.
К примеру, мы хотим в рамках нашей модели проверить, что воздействие, рост фертильности, вызванный увеличением доли мужчин, занятых в сельском хозяйстве, одинаков по силе с ростом фертильности при росте католического населения, то есть я хочу проверить гипотезу о том, что разница этих двух коэффициентов равна нулю. Как это сделать?
Гипотеза не отвергается
Второй способ — подобрать коэффициенты. Т.е. отнять и прибавить коэффициент при втором регрессоре. Программа автоматом рассчитает вероятность того, что разница коэффициентов равна нулю
Гипотеза не отвергается
Подключаем нужные пакеты
library(memisc)
## Loading required package: lattice
## Loading required package: MASS
##
## Attaching package: 'memisc'
## The following object is masked from 'package:ggplot2':
##
## syms
## The following objects are masked from 'package:stats':
##
## contr.sum, contr.treatment, contrasts
## The following object is masked from 'package:base':
##
## as.array
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ tibble 2.1.3 ✔ purrr 0.3.2
## ✔ tidyr 0.8.2 ✔ dplyr 0.8.3
## ✔ readr 1.3.1 ✔ stringr 1.4.0
## ✔ tibble 2.1.3 ✔ forcats 0.4.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%@%() masks memisc::%@%()
## ✖ dplyr::collect() masks memisc::collect()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::recode() masks memisc::recode()
## ✖ dplyr::rename() masks memisc::rename()
## ✖ dplyr::select() masks MASS::select()
## ✖ dplyr::syms() masks memisc::syms(), ggplot2::syms()
library(psych)
##
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(sjPlot)
library(sgof)
## Loading required package: poibin
library(foreign)
library(car)
## Loading required package: carData
##
## Attaching package: 'car'
## The following object is masked from 'package:psych':
##
## logit
## The following object is masked from 'package:dplyr':
##
## recode
## The following object is masked from 'package:purrr':
##
## some
## The following object is masked from 'package:memisc':
##
## recode
library(hexbin)
library(rlms)
100 величин, распределение нормально. Любое распределение генерится в R функцией, которая начинается с r и далее название распределения.
z <- rnorm(100, mean = 5, sd = 3)
z[56]
## [1] 4.414769
z[2:9]
## [1] 6.986693 3.969756 1.714242 3.768766 3.898127 8.314513 6.326293 8.544909
qplot(z)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Все функции плотности начинаются с буквы d (от density) в R.
x <- seq(-10, 15, by = 0.5)
y <- dnorm(x, mean = 5, sd = 3)
qplot(x, y, geom = "line")
Все функции плотности начинаются с буквы p (от probability) в R.
pnorm(3, mean = 5, sd = 3)
## [1] 0.2524925
если я хочу найти вероятность того, что z лежит в диапазоне от 4 до 9, то это есть с точки зрения здравого смысла вероятность того, что z меньше 9 минус вероятность того, что z меньше 4
pnorm(9, mean = 5, sd = 3) - pnorm(4, mean = 5, sd = 3)
## [1] 0.5393474
Все функции плотности начинаются с буквы q (от quantile) в R.
qnorm(0.7, mean = 5, sd = 3)
## [1] 6.573202
Квантиль — это на самом деле, обратная функция к функции распределения. То есть, если, например, я хочу найти такое число а, чтобы вероятность того, что z меньше а была равна, скажем, 0.7. И вот надо найти такое а. То соответственно, это можно найти с помощью квантильной функции
Есть разные популярные распределения. chisq — хи квадрат. t — стьюдента. f — Эф распределение
Множественная регрессия, проверка гипиотез
h <- swiss
glimpse(h)
## Observations: 47
## Variables: 6
## $ Fertility <dbl> 80.2, 83.1, 92.5, 85.8, 76.9, 76.1, 83.8, 92.4,…
## $ Agriculture <dbl> 17.0, 45.1, 39.7, 36.5, 43.5, 35.3, 70.2, 67.8,…
## $ Examination <int> 15, 6, 5, 12, 17, 9, 16, 14, 12, 16, 14, 21, 14…
## $ Education <int> 12, 9, 5, 7, 15, 7, 7, 8, 7, 13, 6, 12, 7, 12, …
## $ Catholic <dbl> 9.96, 84.84, 93.40, 33.77, 5.16, 90.57, 92.85, …
## $ Infant.Mortality <dbl> 22.2, 22.2, 20.2, 20.3, 20.6, 26.6, 23.6, 24.9,…
Мы оценим линейную модель регрессии. Будем предполагать, что фертильность, Fertility, зависит от доли католического населения в данном кантоне, от показателя, насколько это регион сельскохозяйственный, и, скажем, от Examination.
model <- lm(data = h, Fertility ~ Catholic + Agriculture + Examination)
Посмотрим отчет по этой моделе
summary(model)
##
## Call:
## lm(formula = Fertility ~ Catholic + Agriculture + Examination,
## data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.746 -5.724 1.248 6.324 18.849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.86803 8.63691 10.521 1.81e-13 ***
## Catholic 0.04240 0.04148 1.022 0.312347
## Agriculture -0.09516 0.08587 -1.108 0.273930
## Examination -1.07035 0.27316 -3.918 0.000315 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.616 on 43 degrees of freedom
## Multiple R-squared: 0.4461, Adjusted R-squared: 0.4074
## F-statistic: 11.54 on 3 and 43 DF, p-value: 1.113e-05
Отсюда можем сказать, что гипотеза о том, что \(\beta_1\) равно нулю, отвергается; гипотеза о том, что \(\beta_2\) равно нулю, не отвергается; гипотеза о том, что \(\beta_3\) равно нулю, не отвергается; гипотеза о том, что \(\beta_4\) равно нулю, отвергается.
Только значения коэффициентов, И также можем легко получить доверительные интервалы:
coeftest(model)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.868026 8.636909 10.5209 1.807e-13 ***
## Catholic 0.042401 0.041476 1.0223 0.3123467
## Agriculture -0.095159 0.085867 -1.1082 0.2739301
## Examination -1.070347 0.273160 -3.9184 0.0003147 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
confint(model)
## 2.5 % 97.5 %
## (Intercept) 73.45003853 108.28601395
## Catholic -0.04124219 0.12604494
## Agriculture -0.26832536 0.07800811
## Examination -1.62122538 -0.51946784
проверим линейную гипотезу о том, что коэффициент зависимости при доле католического населения и при доле населения, занятого в сельском хозяйстве, одинаковые.
Воспользуемся хитрым способом, в котором мы складываем две переменные — способ с построением вспомогательной регрессии — чтобы проверить гипотезу.
Значок I означает инструкцию для R, что надо трактовать Catholic + Agriculture — «плюс» в прямом смысле.
model_aux <- lm(data = h,
Fertility ~ Catholic + I(Catholic + Agriculture) + Examination)
Выведем отчёт о модели
summary(model_aux)
##
## Call:
## lm(formula = Fertility ~ Catholic + I(Catholic + Agriculture) +
## Examination, data = h)
##
## Residuals:
## Min 1Q Median 3Q Max
## -26.746 -5.724 1.248 6.324 18.849
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.86803 8.63691 10.521 1.81e-13 ***
## Catholic 0.13756 0.09585 1.435 0.158483
## I(Catholic + Agriculture) -0.09516 0.08587 -1.108 0.273930
## Examination -1.07035 0.27316 -3.918 0.000315 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.616 on 43 degrees of freedom
## Multiple R-squared: 0.4461, Adjusted R-squared: 0.4074
## F-statistic: 11.54 on 3 and 43 DF, p-value: 1.113e-05
Коэффициент при Catholic незначим, потому что P-уровень равен 0.158483. Это говорит о том, что гипотеза о том, что коэффициенты при Catholic и Agriculture равны, не отвергается.
Таким образом, мы смогли проверить гипотезу о том, что два коэффициента истинных, неизвестных нам, равны, и эта гипотеза в нашем случае не отвергается.
Один из способов почувствовать существенный коэффициент или нет, — это посчитать стандартизированные коэффициенты \(\widehat\beta\), то есть привести все объясняющие переменные и объясняемую переменную к неким универсальным единицам измерения, чтобы они были сравнимы, а именно: вычесть из каждой переменной её среднее и поделить на оценённое стандартное отклонение.
Шаг 1 — Стандартизация коэффициентов
На шаге один мы преобразуем каждую переменную, масштабируем каждую переменную. Значит, создадим набор h_st, где мы изменим каждую переменную, функция называется mutate_each, в наборе данных h, а формула, по которой мы будем, функция, по которой мы будем менять каждую переменную, называется scale, она осуществляет как раз масштабирование, вычитание среднего и деления на стандартную ошибку.
h_st <- mutate_each(h, "scale")
glimpse(h_st)
## Observations: 47
## Variables: 6
## $ Fertility <dbl[,1]> <matrix[25 x 1]>
## $ Agriculture <dbl[,1]> <matrix[25 x 1]>
## $ Examination <dbl[,1]> <matrix[25 x 1]>
## $ Education <dbl[,1]> <matrix[25 x 1]>
## $ Catholic <dbl[,1]> <matrix[25 x 1]>
## $ Infant.Mortality <dbl[,1]> <matrix[25 x 1]>
Шаг 2 — построение модели
model_st <- lm(data = h_st, Fertility ~ Catholic + Agriculture + Examination)
summary(model_st)
##
## Call:
## lm(formula = Fertility ~ Catholic + Agriculture + Examination,
## data = h_st)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.14113 -0.45821 0.09989 0.50623 1.50891
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.729e-16 1.123e-01 0.000 1.000000
## Catholic 1.416e-01 1.385e-01 1.022 0.312347
## Agriculture -1.730e-01 1.561e-01 -1.108 0.273930
## Examination -6.836e-01 1.745e-01 -3.918 0.000315 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7698 on 43 degrees of freedom
## Multiple R-squared: 0.4461, Adjusted R-squared: 0.4074
## F-statistic: 11.54 on 3 and 43 DF, p-value: 1.113e-05
Шаг 3 — визуализация коэффициентов.
Шаг 4 — принятие решения
Сейчас мы на искусственных данных проиллюстрируем идею того, что нельзя просто так построить регрессию и отвечать на вопрос, а какие коэффициенты у меня значимы.
Мы сочиним в нашем искусственном эксперименте совершенно несвязанный игрек, который никак не зависит от якобы объясняющих переменных. У нас будет 40 якобы объясняющих переменных, одна якобы зависимая, хотя на самом деле независимая, и мы будем строить, оценивать модель линейной регрессии.
Генерим данные
set.seed(42)
d <- matrix(rnorm(100*41, mean = 0, sd = 1), nrow = 100)
df <- data.frame(d)
glimpse(df)
## Observations: 100
## Variables: 41
## $ X1 <dbl> 1.37095845, -0.56469817, 0.36312841, 0.63286260, 0.40426832,…
## $ X2 <dbl> 1.200965376, 1.044751087, -1.003208647, 1.848481902, -0.6667…
## $ X3 <dbl> -2.00092924, 0.33377720, 1.17132513, 2.05953924, -1.37686160…
## $ X4 <dbl> -0.0046207678, 0.7602421677, 0.0389909129, 0.7350721416, -0.…
## $ X5 <dbl> 1.33491259, -0.86927176, 0.05548695, 0.04906691, -0.57835573…
## $ X6 <dbl> 1.029140719, 0.914774868, -0.002456267, 0.136009552, -0.7201…
## $ X7 <dbl> -0.2484829, 0.4223204, 0.9876533, 0.8355682, -0.6605219, 1.5…
## $ X8 <dbl> 0.29469236, 0.39274127, -1.00084371, -0.32572712, -1.0083488…
## $ X9 <dbl> 0.68880775, 0.72508302, 0.21738021, -0.20165673, -1.36568986…
## $ X10 <dbl> 0.94192422, -0.24861404, 0.09647886, -0.43393094, 2.17866787…
## $ X11 <dbl> 2.32505849, 0.52412218, 0.97073342, 0.37697340, -0.99593340,…
## $ X12 <dbl> -0.6502844, -1.0031832, -0.5351139, -0.1104146, 0.6004302, 0…
## $ X13 <dbl> -0.74651645, 0.03660612, 0.32330962, 0.37967603, 0.87655650,…
## $ X14 <dbl> -1.18554694, -0.65838929, 1.08950795, 0.50878594, -0.1359066…
## $ X15 <dbl> 0.87729465, -1.77337144, -0.04568732, -0.39487225, -0.128056…
## $ X16 <dbl> -0.60138300, -0.13581614, -0.98727278, 0.83192501, -0.795059…
## $ X17 <dbl> 0.02931665, 1.93431434, 0.83751088, -0.18534633, 0.53332918,…
## $ X18 <dbl> 1.44820167, -0.08585566, -1.57376980, 1.29247616, 0.08446776…
## $ X19 <dbl> -1.2213334, -0.4529860, -0.7023199, -0.6743804, -2.3030232, …
## $ X20 <dbl> 0.09782704, 1.46502103, 0.99043585, 1.41438675, -1.33589777,…
## $ X21 <dbl> 0.25057807, -0.27792405, -1.72473573, -2.00670494, -1.291808…
## $ X22 <dbl> -0.58693858, 0.92697573, -0.06540585, -0.93711176, -0.633859…
## $ X23 <dbl> 0.19033984, -0.07173877, -0.00285171, -1.10821896, 0.9351917…
## $ X24 <dbl> 0.75806083, 1.47961434, -0.61956241, 0.25968723, -1.08652092…
## $ X25 <dbl> 1.32750456, -0.60083583, 0.05650686, -0.53107629, -0.0808987…
## $ X26 <dbl> 0.617336710, -0.004541141, -0.091256360, 0.399959375, 0.5889…
## $ X27 <dbl> -0.2418159, 0.5721363, 0.7057918, 0.2048996, 0.5521034, -0.5…
## $ X28 <dbl> -0.54468934, -0.43453296, -0.14084519, 1.75814127, -0.830141…
## $ X29 <dbl> -0.8495311, 0.1151387, -1.3313372, 2.3744583, -2.4958364, 2.…
## $ X30 <dbl> -0.63258679, 1.06719020, 0.33739592, -0.95364566, 0.25117357…
## $ X31 <dbl> -0.68566166, -0.79271447, -0.40700415, -1.14867061, 1.115760…
## $ X32 <dbl> 1.37843190, -2.95313731, 0.95583656, 1.06003423, 1.35096057,…
## $ X33 <dbl> -2.0645260, 0.6772492, 0.8363831, 0.1222439, -0.7882496, 0.6…
## $ X34 <dbl> -0.9762936, -1.0083745, -0.2563995, 0.8828167, 0.3760271, 0.…
## $ X35 <dbl> -0.1062227, -1.2870865, -0.9080599, -0.2839357, 0.2665925, 1…
## $ X36 <dbl> 0.29043335, -1.30379423, -0.11410060, -0.09296209, -0.122560…
## $ X37 <dbl> -0.18434938, -0.26304466, -0.52564624, -0.67522750, -0.15678…
## $ X38 <dbl> -0.59179607, -0.88500175, -0.52743656, 0.03103741, -0.946515…
## $ X39 <dbl> -0.007239089, -0.462310513, -0.724387536, -0.639016920, 0.44…
## $ X40 <dbl> -0.4503987, -1.3772278, 1.0485915, 0.3055562, -0.9187486, -0…
## $ X41 <dbl> -0.1418087, -0.8138981, -0.3255406, 0.3781574, -1.9944854, -…
Будем объяснять первую переменную, всеми оставшимися в наборе данных переменными. Для этого есть удобное сокращение в виде точки .
model_pusto <- lm(data = df, X1~.)
summary(model_pusto)
##
## Call:
## lm(formula = X1 ~ ., data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.7476 -0.6018 0.1303 0.5730 1.5573
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001317 0.143255 0.009 0.9927
## X2 0.142212 0.161744 0.879 0.3828
## X3 -0.314385 0.136043 -2.311 0.0244 *
## X4 0.031564 0.166652 0.189 0.8504
## X5 0.178323 0.150741 1.183 0.2416
## X6 0.111487 0.137823 0.809 0.4218
## X7 -0.108252 0.137431 -0.788 0.4340
## X8 0.145668 0.149928 0.972 0.3352
## X9 -0.036021 0.143237 -0.251 0.8023
## X10 0.110686 0.121645 0.910 0.3666
## X11 -0.038508 0.144271 -0.267 0.7905
## X12 0.034769 0.145851 0.238 0.8124
## X13 0.049680 0.152415 0.326 0.7456
## X14 -0.158988 0.142948 -1.112 0.2706
## X15 -0.014374 0.148596 -0.097 0.9233
## X16 0.248442 0.172522 1.440 0.1551
## X17 -0.061389 0.144079 -0.426 0.6716
## X18 -0.047658 0.138050 -0.345 0.7312
## X19 0.235795 0.141078 1.671 0.0999 .
## X20 -0.034408 0.148917 -0.231 0.8181
## X21 0.001259 0.131329 0.010 0.9924
## X22 -0.180049 0.154533 -1.165 0.2487
## X23 -0.072289 0.132548 -0.545 0.5875
## X24 -0.032262 0.190246 -0.170 0.8659
## X25 0.052217 0.148228 0.352 0.7259
## X26 -0.206708 0.139283 -1.484 0.1431
## X27 -0.082317 0.121478 -0.678 0.5007
## X28 -0.096665 0.164121 -0.589 0.5581
## X29 0.202307 0.133146 1.519 0.1340
## X30 -0.047584 0.129551 -0.367 0.7147
## X31 0.055963 0.162809 0.344 0.7323
## X32 0.307859 0.158290 1.945 0.0566 .
## X33 0.044117 0.150258 0.294 0.7701
## X34 -0.093703 0.135581 -0.691 0.4922
## X35 -0.344475 0.168765 -2.041 0.0457 *
## X36 -0.065275 0.135642 -0.481 0.6321
## X37 0.063727 0.143506 0.444 0.6586
## X38 0.132398 0.161363 0.820 0.4152
## X39 -0.055207 0.158315 -0.349 0.7285
## X40 -0.059821 0.174217 -0.343 0.7325
## X41 -0.085292 0.142895 -0.597 0.5529
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.123 on 59 degrees of freedom
## Multiple R-squared: 0.3074, Adjusted R-squared: -0.1622
## F-statistic: 0.6546 on 40 and 59 DF, p-value: 0.9211
У нас получилось, что переменные X3 и X35 ылияют на первую переменную при 5% уровне значимости, а X19 и X32 влияют на 10% уровне значимости.
С чем это связано? Это связано с тем, что когда мы фиксируем уровень значимости, мы соглашаемся на некоторую вероятность ошибиться.
Соответственно, когда мы зафиксировали уровень значимости 10%, это означает, что с вероятностью 10% мы в случае на самом деле отсутствия зависимости якобы её обнаружим.
Соответственно, из этого искусственного эксперимента нужно сделать простой вывод, что стратегия «я оценю модель с большим количеством объясняющих переменных и выпишу из них значимые, и скажу, что от них игрек зависит», — неправильная, потому что в силу того, что есть для каждого регрессора десятипроцентный или пятипроцентный шанс сделать ошибку, при большом количестве регрессоров кто-то якобы значимым будет, даже если на самом деле никакой зависимости нет.
Есть несколько моделей
model <- lm(data = h, Fertility ~ Catholic + Agriculture + Examination)
model2 <- lm(data = h, Fertility ~ Catholic + Agriculture)
Как нам лучше их сравнить. Сделаем табличку, сравниваем при помощи функции mtable() из пакета memisc
compare_12 <- mtable(model, model2)
beta <- c(500, 30)